#ifdef MPE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <athread.h>
#include <mpi.h>
#include "force_lj_sw64.h"
#include "swcache.h"
//#include "types.h"
//#include <mpi.h>

extern SLAVE_FUN(compute_original_cpe)(compute_param *);
extern SLAVE_FUN(compute_fullneigh_cpe_a2s)(compute_param *);
extern SLAVE_FUN(compute_fullneigh_cpe)(compute_param *);
extern SLAVE_FUN(compute_cpe_a2s)(compute_param *);

void compute_fullneigh_c(compute_param* param)
{
  int proc;
  MPI_Comm_rank(MPI_COMM_WORLD,&proc);
  
  int    i,k;
  int    nlocal      = param->nlocal;
  int    nall        = param->nall;
  int    *type       = param->type;
  double *x          = param->x;
  double *f          = param->f;
  double eng_vdwl    = param->eng_vdwl;
  double virial      = param->virial;
  int    *neighbors  = param->neighbors;
  int    *numneigh   = param->numneigh;
  int    maxneighs   = param->maxneighs;
  double *sigma6     = param->sigma6;
  double *epsilon    = param->epsilon;
  int     EVFLAG     = param->EVFLAG;
  double *cutforcesq = param->cutforcesq;
  int    ntypes      = param->ntypes;
  double t_eng_vdwl = 0;
  double t_virial = 0;
  param->proc = proc;

  //if(proc == 0)
  //  printf("maxneighs = %d\n",maxneighs);//100

//#define SUNWAY
#ifdef SUNWAY
  if(athread_idle() == 0)
    athread_init();

  atom_in_t *atom_in = malloc(sizeof(atom_in_t)*(nall + 32));
  param->atom_in = (void*)((long)atom_in | 255) + 1;

  athread_spawn(compute_fullneigh_cpe_a2s,param);
  athread_join();

  //atom_in_t* atom_test = param->atom_in;

  //for(i = 0;i < nlocal;i++)
  //{
  //  if(atom_test[i].x[0] - x[i*3+0] > 0.000001 || \ 
  //     atom_test[i].x[1] - x[i*3+1] > 0.000001 || \
  //     atom_test[i].x[2] - x[i*3+2] > 0.000001 || \
  //     atom_test[i].type != type[i])
  //     printf("fsfsafdsfdsafsdafsadfasd\n");
  //}
  
  athread_spawn(compute_fullneigh_cpe,param);
  athread_join();

#else

  for(i = 0; i < nlocal; i++) {
    f[i * PAD + 0] = 0.0;
    f[i * PAD + 1] = 0.0;
    f[i * PAD + 2] = 0.0;
  }

  for(i = 0; i < nlocal; i++) {
    const int* const neighs = &neighbors[i * maxneighs];
    const int numneighs = numneigh[i];
    const double xtmp = x[i * PAD + 0];
    const double ytmp = x[i * PAD + 1];
    const double ztmp = x[i * PAD + 2];
    const int type_i = type[i];
    double fix = 0;
    double fiy = 0;
    double fiz = 0;

    for(k = 0; k < numneighs; k++) {
      const int j = neighs[k];
      const double delx = xtmp - x[j * PAD + 0];
      const double dely = ytmp - x[j * PAD + 1];
      const double delz = ztmp - x[j * PAD + 2];
      const int type_j = type[j];
      const double rsq = delx * delx + dely * dely + delz * delz;

      int type_ij = type_i*ntypes+type_j;
      
      if(rsq < cutforcesq[type_ij]) {
        const double sr2 = 1.0 / rsq;
        const double sr6 = sr2 * sr2 * sr2 * sigma6[type_ij];
        const double force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[type_ij];
        fix += delx * force;
        fiy += dely * force;
        fiz += delz * force;

        if(EVFLAG) {
          t_eng_vdwl += sr6 * (sr6 - 1.0) * epsilon[type_ij];
          t_virial += (delx * delx + dely * dely + delz * delz) * force;
        }
      }
      
    }

    f[i * PAD + 0] += fix;
    f[i * PAD + 1] += fiy;
    f[i * PAD + 2] += fiz;

  }

  t_eng_vdwl *= 4.0;
  t_virial *= 0.5;

  eng_vdwl += t_eng_vdwl;
  virial += t_virial;

  param->eng_vdwl = eng_vdwl;
  param->virial   = virial;

  //if(proc == 0)
  //  printf("eng_vdwl = %lf,virial = %lf\n",eng_vdwl,virial);

#endif
}

swcache_lock_t *locks = NULL;
void compute_original_c(compute_param* param)
{
  int proc;
  MPI_Comm_rank(MPI_COMM_WORLD,&proc);
  
  int    i,k;
  int    nlocal      = param->nlocal;
  int    nall        = param->nall;
  int    *type       = param->type;
  double *x          = param->x;
  double *f          = param->f;
  double eng_vdwl    = param->eng_vdwl;
  double virial      = param->virial;
  int    *neighbors  = param->neighbors;
  int    *numneigh   = param->numneigh;
  int    maxneighs   = param->maxneighs;
  double *sigma6     = param->sigma6;
  double *epsilon    = param->epsilon;
  int     EVFLAG     = param->EVFLAG;
  double *cutforcesq = param->cutforcesq;
  int    ntypes      = param->ntypes;

#define SUNWAY
#ifdef SUNWAY
  //if(proc == 0)
  //  printf("----------CPE--------\n");
#ifndef W_CACHE
  for(i = 0; i < nall; i++) {
    f[i * PAD + 0] = 0.0;
    f[i * PAD + 1] = 0.0;
    f[i * PAD + 2] = 0.0;
  }
#endif

  if(athread_idle() == 0)
    athread_init();
  double start1,start2,start3=0;
  double stop1,stop2,stop3=0;
  atom_in_t *atom_in = malloc(sizeof(atom_in_t)*(nall + 32));
  param->atom_in = (void*)((long)atom_in | 255) + 1;
  
  force_t *force_tmp = malloc(sizeof(force_t)*(nall + 32));
  //memset(force_tmp, 0, sizeof(force_t) * (nall + 32));
  param->force_tmp = (void*)((long)force_tmp | 255) + 1;

  start1 = MPI_Wtime();
  athread_spawn(compute_fullneigh_cpe_a2s,param);
  athread_join();
  stop1 = MPI_Wtime();
  
  
  if (locks == NULL) {
    locks = swcache_u_prepare_locks(nall * 4,5);
  }
  param->locks = locks;
  
  start2 = MPI_Wtime();
  athread_spawn(compute_original_cpe,param);
  athread_join();
  stop2 = MPI_Wtime();

#define W_CACHE
#ifdef W_CACHE
  start3 = MPI_Wtime();
  athread_spawn(compute_cpe_a2s,param);
  athread_join();
  stop3 = MPI_Wtime();
#endif

  if(proc == 0)
    printf("a2s compute s2a time = %lf,%lf,%lf\n",stop1 - start1,stop2 - start2,stop3 - start3);

#else
  for(i = 0; i < nall; i++) {
    f[i * PAD + 0] = 0.0;
    f[i * PAD + 1] = 0.0;
    f[i * PAD + 2] = 0.0;
  }

  for(i = 0; i < nlocal; i++) {
    int* neighs = &neighbors[i * maxneighs];
    int numneigh_c = numneigh[i];
    double xtmp = x[i * PAD + 0];
    double ytmp = x[i * PAD + 1];
    double ztmp = x[i * PAD + 2];
    int type_i = type[i];

    for(k = 0; k < numneigh_c; k++) {
      int j = neighs[k];
      double delx = xtmp - x[j * PAD + 0];
      double dely = ytmp - x[j * PAD + 1];
      double delz = ztmp - x[j * PAD + 2];
      int type_j = type[j];
      double rsq = delx * delx + dely * dely + delz * delz;

      int type_ij = type_i*ntypes+type_j;
  
      if(rsq < cutforcesq[type_ij]) {
        double sr2 = 1.0 / rsq;
        double sr6 = sr2 * sr2 * sr2 * sigma6[type_ij];
        double force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[type_ij];
       
        f[i * PAD + 0] += delx * force;
        f[i * PAD + 1] += dely * force;
        f[i * PAD + 2] += delz * force;
        f[j * PAD + 0] -= delx * force;
        f[j * PAD + 1] -= dely * force;
        f[j * PAD + 2] -= delz * force;

        if(EVFLAG) {
          eng_vdwl += (4.0 * sr6 * (sr6 - 1.0)) * epsilon[type_ij];
          virial += (delx * delx + dely * dely + delz * delz) * force;
        }
      }
    }
  }
  param->eng_vdwl = eng_vdwl;
  param->virial   = virial;
#endif
}

#endif

#ifdef CPE
#include <slave.h>
#include <simd.h>
#include "force_lj_sw64.h"
#include <dma.h>
#include "/home/export/online1/systest/swsduhpc/dxh/workspace/miniFE-ref/src/dma_macros.h"
//#include "/home/export/online1/systest/swsduhpc/dxh/workspace/miniFE-ref/src/swcache.h"
#include "swcache.h"


#define CACHE_HBIT 6
#define CACHE_SBIT 3
#define CACHE_LINESIZE (1 << CACHE_SBIT)
#define CACHE_LINECNT (1 << CACHE_HBIT)
#define CACHE_MMASK (CACHE_LINESIZE - 1)
#define CACHE_LMASK (CACHE_LINECNT - 1)


void reg_reduce_inplace_doublev4(doublev4 *arr, int len)
{
  int i, j, n;
  doublev4 tmp;
  for (i = 1; i < 8; i += i)
  {
    if ((_ROW & i) == i)
    {
      for (j = 0; j < len; j ++)
        asm("putc %0, %1": : "r"(arr[j]), "r"(_ROW ^ i));
    }
    if ((_ROW & i) == 0)
    {
      for (j = 0; j < len; j ++)
      {
        asm("getc %0\n\t"
            "vaddd %0, %1, %1\n\t"
            : "=r"(tmp), "+r"(arr[j]));
      }
    }
    athread_syn(COL_SCOPE, 0xff);
  }

  athread_syn(ARRAY_SCOPE, 0xffff);
  if (_ROW == 0)
  {
    for (i = 1; i < 8; i += i)
    {
      if ((_COL & i) == i)
      {
        for (j = 0; j < len; j ++)
          asm("putr %0, %1": : "r"(arr[j]), "r"(_COL ^ i));
      }
      if ((_COL & i) == 0)
      {
        for (j = 0; j < len; j ++)
        {
          asm("getr %0\n\t" 
              "vaddd %0, %1, %1\n\t"
              : "=r"(tmp), "+r"(arr[j]));
        }
      }
    }
    athread_syn(ROW_SCOPE, 0xff);
  }
}

void compute_cpe_a2s(compute_param *param)
{
  dma_init();
  compute_param lp;
  pe_get(param,&lp,sizeof(compute_param));
  dma_syn();

  int nlocal = lp.nlocal;
  int nall   = lp.nall;

  double f[ISTEP*3];
  force_t force_tmp[ISTEP];
  int ist,ied,isz,i,j,k;

  for(ist = _MYID*ISTEP;ist < nall;ist += ISTEP*64)
  {
    ied = ist + ISTEP;
    if(ied > nall)
      ied = nall;
    isz = ied - ist;
    pe_get(lp.force_tmp + ist, force_tmp, sizeof(force_t)*isz);
    dma_syn();
    for(i = 0;i < isz;i++)
    {
      f[PAD*i + 0] = force_tmp[i].f[0];
      f[PAD*i + 1] = force_tmp[i].f[1];
      f[PAD*i + 2] = force_tmp[i].f[2];
    }
    pe_put(lp.f+ist*PAD,f,sizeof(double)*isz*3);
    dma_syn();
  }
}

void compute_fullneigh_cpe_a2s(compute_param *param)
{
  dma_init();
  compute_param lp;
  pe_get(param,&lp,sizeof(compute_param));
  dma_syn();

  int nlocal = lp.nlocal;
  int nall   = lp.nall;

  double x[ISTEP*3];
  int type[ISTEP];
  atom_in_t atom_in[ISTEP];
  int ist,ied,isz,i,j,k;
  force_t force_tmp[ISTEP];

  for(ist = _MYID*ISTEP;ist < nall;ist += ISTEP*64)
  {
    ied = ist + ISTEP;
    if(ied > nall)
      ied = nall;
    isz = ied - ist;
    pe_get(lp.x + ist*3, x, sizeof(double)*isz*3);
    pe_get(lp.type + ist, type, sizeof(int)*isz);
    dma_syn();
    for(i = 0;i < isz;i++)
    {
      atom_in[i].x[0] =  x[i*3+0];
      atom_in[i].x[1] =  x[i*3+1];
      atom_in[i].x[2] =  x[i*3+2];
      atom_in[i].type = type[i];
      force_tmp[i].f[0] = 0.0;
      force_tmp[i].f[1] = 0.0;
      force_tmp[i].f[2] = 0.0;
    }
    pe_put(lp.atom_in+ist,atom_in,sizeof(atom_in_t)*isz);
    pe_put(lp.force_tmp+ist,force_tmp,sizeof(force_t)*isz);
    dma_syn();
  }
}

void compute_fullneigh_cpe(compute_param *param)
{

  dma_init();
  compute_param lp;
  pe_get(param,&lp,sizeof(compute_param));
  dma_syn();

  int i,j,k,ii,jj,kk;
  int    nlocal      = lp.nlocal;
  int    nall        = lp.nall;
  int    *type       = lp.type;
  double *x          = lp.x;
  double *f          = lp.f;
  double eng_vdwl    = lp.eng_vdwl;
  double virial      = lp.virial;
  int    *neighbors  = lp.neighbors;
  int    *numneigh   = lp.numneigh;
  int    maxneighs   = lp.maxneighs;
  int     EVFLAG     = lp.EVFLAG;
  //double *cutforcesq = lp.cutforcesq;
  int    ntypes      = lp.ntypes;
  double t_eng_vdwl = 0;
  double t_virial = 0;
  int    proc = lp.proc;
  atom_in_t* atom_jj = lp.atom_in;

  int ist,ied,isz;
  int jst,jed,jsz;

  doublev4 eng_virial_v4 = 0;
  double *vdwl = &eng_virial_v4;
  double *virl = vdwl + 1;

  double    fi[ISTEP*PAD];
  double    xi[ISTEP*PAD];
  int       typei[ISTEP];
  int       numneighi[ISTEP];
  atom_in_t atom_in[ISTEP];
  int    neighsi[JSTEP];
  double INF = 1.0/0.0;

//#define SWCACHE
#ifndef SWCACHE
  atom_in_t atom_buff[JSTEP];
  atom_in_t atom_cache[CACHE_LINECNT][CACHE_LINESIZE];
  int       tag_cache[CACHE_LINECNT];

  for(i = 0;i < CACHE_LINECNT;i++)
    tag_cache[i] = -1;
#else
#define MINI_CACHE atom_jj, atom_in_t, 5, 3, 7 
  SWCACHE_INIT(MINI_CACHE);
#endif

  double sigma6[ntypes*ntypes];
  double epsilon[ntypes*ntypes];
  double cutforcesq[ntypes*ntypes];

  pe_get(lp.sigma6,sigma6,sizeof(double)*ntypes*ntypes);
  pe_get(lp.epsilon,epsilon,sizeof(double)*ntypes*ntypes);
  pe_get(lp.cutforcesq,cutforcesq,sizeof(double)*ntypes*ntypes);
  dma_syn();


  for(ist = _MYID * ISTEP; ist < nlocal; ist += ISTEP * 64) {
    
    ied = ist + ISTEP;
    if(ied > nlocal)
      ied = nlocal;
    isz = ied - ist;

    pe_get(numneigh + ist, numneighi,sizeof(int)*isz);
    pe_get(lp.atom_in + ist, atom_in,sizeof(atom_in_t)*isz);
    dma_syn();
    
    for(i = ist;i < ied;i++)
    {
      int ioff = i - ist;

      fi[ioff*PAD+0] = 0.0;
      fi[ioff*PAD+1] = 0.0;
      fi[ioff*PAD+2] = 0.0;

      int* neighs = neighbors + i * maxneighs;
      int numneighs = numneighi[ioff];
      double xtmp = atom_in[ioff].x[0];
      double ytmp = atom_in[ioff].x[1];
      double ztmp = atom_in[ioff].x[2];
      
      int type_i = atom_in[ioff].type;

      double fix = 0;
      double fiy = 0;
      double fiz = 0;

      for(jst = 0; jst < numneighs; jst += JSTEP )
      {
        jed = jst + JSTEP;
        if(jed > numneighs)
          jed = numneighs;
        jsz = jed - jst;
        pe_get(neighs + jst,neighsi,sizeof(int)*jsz);
        dma_syn();

        for(jj = 0;jj < jsz; jj++)
        {
#ifndef SWCACHE
          j = neighsi[jj];
          int cache_tag = j >> CACHE_SBIT;
          int cache_line = j >> CACHE_SBIT & CACHE_LMASK;
          int cache_head = j &~ CACHE_MMASK;
          int cache_off = j & CACHE_MMASK;

          if(tag_cache[cache_line] != cache_tag)
          {
            pe_get(lp.atom_in + cache_head,atom_cache[cache_line],sizeof(atom_in_t)*CACHE_LINESIZE);
            dma_syn();
            tag_cache[cache_line] = cache_tag;
          }
          atom_buff[jj] = atom_cache[cache_line][cache_off];
          
          double delx = xtmp - atom_buff[jj].x[0];
          double dely = ytmp - atom_buff[jj].x[1];
          double delz = ztmp - atom_buff[jj].x[2];
          int type_j = atom_buff[jj].type;
#else
          atom_in_t *atom_data;
          int idx = jj;
          SWCACHE_TADDR(atom_data,MINI_CACHE,neighsi[jj]);

          double delx = xtmp - atom_data->x[0];
          double dely = ytmp - atom_data->x[1];
          double delz = ztmp - atom_data->x[2];
          int type_j = atom_data->type;
#endif
          double rsq = delx * delx + dely * dely + delz * delz;
          int type_ij = type_i*ntypes+type_j;
  
          if(rsq < cutforcesq[type_ij]) {
            const double sr2 = 1.0 / rsq;
            const double sr6 = sr2 * sr2 * sr2 * sigma6[type_ij];
            const double force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[type_ij];
            fix += delx * force;
            fiy += dely * force;
            fiz += delz * force;
            
            if(EVFLAG) {
              t_eng_vdwl += sr6 * (sr6 - 1.0) * epsilon[type_ij];
              t_virial += (delx * delx + dely * dely + delz * delz) * force;
            }
          }
        }
      
      }

      fi[ioff * PAD + 0] += fix;
      fi[ioff * PAD + 1] += fiy;
      fi[ioff * PAD + 2] += fiz;
    }
    pe_put(f+ist*3,fi,sizeof(double)*isz*3);
    dma_syn();
  }
    
  vdwl[0] = t_eng_vdwl;
  vdwl[1] = t_virial;

  reg_reduce_inplace_doublev4(&eng_virial_v4,1);

  if(_MYID == 0)
  {
    vdwl[0] *= 4.0;
    vdwl[1] *= 0.5;

    eng_vdwl += vdwl[0];
    virial   += vdwl[1];

    vdwl[0] = eng_vdwl;
    vdwl[1] = virial;

    pe_put(&(param->eng_vdwl),&eng_virial_v4,sizeof(double) * 2);
    dma_syn();
  }

  //if(_MYID == 0 && lp.proc == 0)
  //  printf("eng_vdwl = %lf,virial = %lf\n",eng_vdwl,virial);

}


void compute_original_cpe(compute_param *param)
{
  dma_init();
  compute_param lp;
  pe_get(param,&lp,sizeof(compute_param));
  dma_syn();

  int i,j,k,ii,jj,kk;
  int    nlocal      = lp.nlocal;
  int    nall        = lp.nall;
  int    *type       = lp.type;
  double *x          = lp.x;
  double *f          = lp.f;
  double eng_vdwl    = lp.eng_vdwl;
  double virial      = lp.virial;
  int    *neighbors  = lp.neighbors;
  int    *numneigh   = lp.numneigh;
  int    maxneighs   = lp.maxneighs;
  int     EVFLAG     = lp.EVFLAG;
  int    ntypes      = lp.ntypes;
  double t_eng_vdwl = 0;
  double t_virial = 0;
  //param->proc = proc;
  atom_in_t* atom_jj = lp.atom_in;
  force_t *force_tmp = lp.force_tmp;
  swcache_lock_t *locks = lp.locks;

  int ist,ied,isz;
  int jst,jed,jsz;

  doublev4 eng_virial_v4 = 0;
  double *vdwl = &eng_virial_v4;
  double *virl = vdwl + 1;

  double    fi[ISTEP*PAD];
  //double    fj[ISTEP*PAD];
  double    xi[ISTEP*PAD];
  int       typei[ISTEP];
  int       numneighi[ISTEP];
  atom_in_t atom_in[ISTEP];
  int    neighsi[JSTEP];
  double INF = 1.0/0.0;
  //int       tag_cache_w[CACHE_LINECNT];
  
  //for(i = 0;i < CACHE_LINECNT;i++)
  //tag_cache_w[i] = -1;

#define SWCACHE
#ifndef SWCACHE
  atom_in_t atom_buff[JSTEP];
  atom_in_t atom_cache[CACHE_LINECNT][CACHE_LINESIZE];
  int       tag_cache[CACHE_LINECNT];

  for(i = 0;i < CACHE_LINECNT;i++)
    tag_cache[i] = -1;
#else
#define MINI_CACHE atom_jj, atom_in_t, 5, 3, 6
  SWCACHE_INIT(MINI_CACHE);
#endif

#define DT_INIT_OP(x) {*x = 0;}
#define DT_MERGE_OP(x,y) {*x += *y;}

#define W_CACHE
#ifdef W_CACHE
#define DT_CACHE force_tmp,force_t,5,3,6,FADDD
#else
#define DT_CACHE f,double,3,5,6,FADDD
#endif

  SWCACHE_INIT_U(DT_CACHE,locks);

  double sigma6[ntypes*ntypes];
  double epsilon[ntypes*ntypes];
  double cutforcesq[ntypes*ntypes];

  pe_get(lp.sigma6,sigma6,sizeof(double)*ntypes*ntypes);
  pe_get(lp.epsilon,epsilon,sizeof(double)*ntypes*ntypes);
  pe_get(lp.cutforcesq,cutforcesq,sizeof(double)*ntypes*ntypes);
  dma_syn();
  int cpeid;
  asm ("rcsr %0, 0\n\t" : "=r"(cpeid));
  for(ist = _MYID * ISTEP; ist < nlocal; ist += ISTEP * 64) {
    
    ied = ist + ISTEP;
    if(ied > nlocal)
      ied = nlocal;
    isz = ied - ist;
    
    pe_get(numneigh + ist, numneighi,sizeof(int)*isz);
    pe_get(lp.atom_in + ist, atom_in,sizeof(atom_in_t)*isz);
    dma_syn();

    for(i = ist;i < ied;i++)
      {
        int ioff = i - ist;

        int* neighs = neighbors + i * maxneighs;
        int numneighs = numneighi[ioff];
        double xtmp = atom_in[ioff].x[0];
        double ytmp = atom_in[ioff].x[1];
        double ztmp = atom_in[ioff].x[2];
       
        int type_i = atom_in[ioff].type;
      
        fi[ioff*PAD+0] = 0.0;
        fi[ioff*PAD+1] = 0.0;
        fi[ioff*PAD+2] = 0.0;
      
        double fix = 0;
        double fiy = 0;
        double fiz = 0;
       
        //int* neighs = &neighbors[i * maxneighs];
        //int numneigh_c = numneigh[i];

        //double xtmp = x[i * PAD + 0];
        //double ytmp = x[i * PAD + 1];
        //double ztmp = x[i * PAD + 2];
        //int type_i = type[i];

        for(jst = 0; jst < numneighs; jst += JSTEP )
          {
         
            jed = jst + JSTEP;
            if(jed > numneighs)
              jed = numneighs;
            jsz = jed - jst;
            pe_get(neighs + jst,neighsi,sizeof(int)*jsz);
            dma_syn();
         
            for(jj = 0;jj < jsz; jj++)
              {
                int j = neighsi[jj];
#ifndef SWCACHE
                int cache_tag = j >> CACHE_SBIT;
                int cache_line = j >> CACHE_SBIT & CACHE_LMASK;
                int cache_head = j &~ CACHE_MMASK;
                int cache_off = j & CACHE_MMASK;

                if(tag_cache[cache_line] != cache_tag)
                  {
                    pe_get(lp.atom_in + cache_head,atom_cache[cache_line],sizeof(atom_in_t)*CACHE_LINESIZE);
                    dma_syn();
                    tag_cache[cache_line] = cache_tag;
                  }
                atom_buff[jj] = atom_cache[cache_line][cache_off];
          
                double delx = xtmp - atom_buff[jj].x[0];
                double dely = ytmp - atom_buff[jj].x[1];
                double delz = ztmp - atom_buff[jj].x[2];
                int type_j = atom_buff[jj].type;
#else
                atom_in_t *atom_data;
                int idx = jj;
                SWCACHE_TADDR(atom_data,MINI_CACHE,j);

                double delx = xtmp - atom_data->x[0];
                double dely = ytmp - atom_data->x[1];
                double delz = ztmp - atom_data->x[2];
                int type_j = atom_data->type;
#endif
         
                //int j = neighs[jj];
                //double delx = xtmp - x[j * PAD + 0];
                //double dely = ytmp - x[j * PAD + 1];
                //double delz = ztmp - x[j * PAD + 2];
                //int type_j = type[j];
          
                double rsq = delx * delx + dely * dely + delz * delz;
                int type_ij = type_i*ntypes+type_j;
                //if (cpeid==4 && type_ij >= 16) printf("%d %p %d %d %d %d\n", atom_data->type, cutforcesq, type_i, type_j, ntypes, i);
                if(rsq < cutforcesq[type_ij]) {
                  double sr2 = 1.0 / rsq;
                  double sr6 = sr2 * sr2 * sr2 * sigma6[type_ij];
                  double force = 48.0 * sr6 * (sr6 - 0.5) * sr2 * epsilon[type_ij];
          
                  fix += delx * force;
                  fiy += dely * force;
                  fiz += delz * force;

                  double fjx = -delx * force;
                  double fjy = -dely * force;
                  double fjz = -delz * force;
                  /* if (cpeid == 0 && i == 6165) { */
                  /*   printf("%d %d\n", type_i, j); */
                  /* } */
#ifdef W_CACHE
                  force_t* force_dataj;
                  force_dataj->f[0] = fjx;
                  force_dataj->f[1] = fjy;
                  force_dataj->f[2] = fjz;
                  SWCACHE_UPDATE(DT_CACHE,(j),(force_dataj));

#else
                  SWCACHE_UPDATE(DT_CACHE,(j*PAD+0),(&fjx));
                  SWCACHE_UPDATE(DT_CACHE,(j*PAD+1),(&fjy));
                  SWCACHE_UPDATE(DT_CACHE,(j*PAD+2),(&fjz));
#endif
                  /* if (cpeid == 0 && i == 6165) { */
                  /*   printf("%d %d\n", type_i, j); */
                  /* } */

                  //f[j * PAD + 0] -= delx * force;
                  //f[j * PAD + 1] -= dely * force;
                  //f[j * PAD + 2] -= delz * force;

                  if(EVFLAG) {
                    t_eng_vdwl += (4.0 * sr6 * (sr6 - 1.0)) * epsilon[type_ij];
                    t_virial += (delx * delx + dely * dely + delz * delz) * force;
                  }
                }
              }//jj
          } 
#ifdef W_CACHE
        force_t* force_datai;
        force_datai->f[0] = fix;
        force_datai->f[1] = fiy;
        force_datai->f[2] = fiz;
        SWCACHE_UPDATE(DT_CACHE,(i),(force_datai));
#else
        SWCACHE_UPDATE(DT_CACHE, ((i)*PAD+0), &fix);
        SWCACHE_UPDATE(DT_CACHE, ((i)*PAD+1), &fiy);
        SWCACHE_UPDATE(DT_CACHE, ((i)*PAD+2), &fiz);
#endif
        /* fi[ioff * PAD + 0] += fix; */
        /* fi[ioff * PAD + 1] += fiy; */
        /* fi[ioff * PAD + 2] += fiz; */
      }
    /* pe_put(f+ist*3,fi,sizeof(double)*isz*3); */
    /* dma_syn(); */
  }
  SWCACHE_FLUSH(DT_CACHE);

  vdwl[0] = t_eng_vdwl;
  vdwl[1] = t_virial;

  reg_reduce_inplace_doublev4(&eng_virial_v4,1);

  if(_MYID == 0)
    {
      eng_vdwl += vdwl[0];
      virial   += vdwl[1];

      vdwl[0] = eng_vdwl;
      vdwl[1] = virial;

      pe_put(&(param->eng_vdwl),&eng_virial_v4,sizeof(double) * 2);
      dma_syn();
    }
}


#endif












